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Abstract 


In two-phase natural circulation loop, the circulating fluid removes the heat fi-om the 
heat source. With natural circulation being the designed heat removal mode, the flow 
may experience various types of thermohydraulic instabilities, such as flow region 
transition, boiling crisis, Ledinegg instability, density wave oscillations. Pressure 
drop oscillations, acoustic oscillations etc. However most common among these 
instabilities is density wave oscillation. So present work has analyzed this instability. 

For this purpose computer code TINFLOS-D (based on drift flux mode) has been 
obtained by modifying the computer code TINFLOS (based on homogeneous 
model). Fhc code TINFLOS-D solves the four conservations equations using 
Linear-stability analytical technique. The result obtained from the code has been 
validated for stability with four single channel two-phase natural circulation loop. 
Two loops HPNCL and Apsara arc in BARC. Code has also been validated with two 
other natural -circulation loops reported in literature (Furuteara’s and Bergles). 
HPNCL has been validated for steady state also. Parametric sensitivity analysis has 
been done for low quality HPNCL for lower threshold of instability. It has been 
found that natural circulation flow rate first increases then decreases with increase in 
system pressure. The instability reduces with increases in system pressure, power 
and drift velocity. The instability also reduces with increase in core inlet resistance 
and core outlet resistance. The instability increases with increase in riser height, 
downcomer diameter and inlet subcooling. 
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Nomenclature 


A Cross sectional area (m^) 

Q kinematic wave velocity (m/s) 

D hydraulic diameter (m) 

f darcy friction factor 

g acceleration due to gravity (m/s^) 

h enthalpy (J/kg) 

hfg latent heat of vaporization (J/kg) 

j volumetric flux density of mixture (m/s) 

j j- volumetric flux density of liquid (m/s) 

/g volumetric flux density of vapour (m/s) 

K, loss coefficient 

L length of the section (m) 
p pressure (bar) 

Ph heated perimeter(m) 
heat flux (W/m^) 

s stability parameter 
t time (s) 

T temperature (K) 

V specific volume (m^/kg) 

Vfi velocity of single phase liquid (m/s) 

Vjg- drift velocity of vapour (m/s) 

Vm mass velocity of centre of mass of the 

mixture (m/s) 
w mass flow rate (kg/s) 

z axial distance (m) 


Greek Symbols 

a void fraction 

T residence time of fluid in single phase 

region (s) 

Fg vapour generation rate (kg/m -s) 

/I non boiling length in core (m) 

A difference 
p density (kg/m^) 

O phase 

Subscripts 

c core 
f liquid 
g vapour 
m mixture 
r relative 
s single phase 
ss steady state 
sub subcooling 



CHAPTER 1 


LITERATURE REVIEW 


1.1 Introduction 

In a natural circulation loop the circulating fluid removes the heat from a heat 
source and transports it to heat sink, the fluid circulation being established due to the 
buoyancy force. The buoyancy force is the driving force. Buoyancy force is result of 
thermally induced density differences in a gravitational field. The heat sink is located 
at the higher elevation than the heat source. Such loops find application in many 
engineering fields such as, gas turbine blade cooling, solar water heater, geothermal 
systems, electrical machine rotor cooling, transformer cooling, nuclear reactor core 
cooling etc. 

Present generation boiling water reactors utilize forced circulation loops for 
heat transfer from core to the steam generator. But this method has the following 
disadvantages: 

1 . Usage of pump is costly 

2. If pump break down occurs then dissipation of fission heat is hampered. This 
results in tremendous accumulation of heat and may cause core melt down in the 
reactor, and 

3. Back up safety measures have to be used. 

Hence next generation of boiling water reactors should have natural circulation of 
primary coolant under normal operational situations. Fig. 1.1 shows the schematic 
diagram of a natural circulation loop. The loop essentially consists of a tubular 
heated section, a riser, a steam drum or a separator, a condenser, a downcomer, 
and upper- lower horizontal sections. The steam collected in separator or steam 
drum goes to condenser, where it is condensed and then sent back to the separator. 
This loop differs from forced circulation loop only by the absence of pump. 
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Fig 1.1 Schematic of Natural Circulation Loop 



In a two-phase natural circulation loop since heat is removed by natural circulation, 
it is inherently safe. This method of cooling enhances passive safety. 

Many proposed advanced designs of nuclear power reactors use 
two-phase natural circulation of primary coolants. However a major problem with 
these loops is occurrence of flow instabilities. Flow instabilities are undesirable in 
boiling, condensing and other two-phase flow systems for several reasons. 

Boure et. al. [1] divided these instabilities into two classes static instability and 
dynamic instability. 

Static instability is one where a small change in one of the independent 
variables leads to a large change in one of the dependent variables. Examples of 
static instability are listed below: 

(a) Flow regime transition: Cyclic flow pattern transitions and flow pattern 
variation occurs. 

(b) Boiling crisis: Bubbly flow has less void but higher Ap than that of annular 
flow characterized by the Ineffective removal of heat from the heated surface leads 
to wall temperature excursion and flow oscillation . 

(c) The Ledinegg instability or flow excursion instability : Flow undergoes 
sudden, large amplitude excursion to a new, stable operating condition, due to the 
internal pressure gradient being less than the external pressure gradient (imposed 
by the pump). 

(d) Geysering: an erratic boiling region exists, where the the surface temperature 
moves between boiling and natural circulation in irregular cycles, due to the 
presence of gas; effect disappears at higher heat fluxes and higher pressures 

Dynamic instability: a flow is subjected to dynamic instability when inertia and other 
feedback effects have an essential part in the process. 



(a) Density wave oscillations: Occurs due to delay and feedback effects between 
flow rate, density and Ap . 

(b) Pressure drop oscillations: These are very low frequency (0.1 Hz) oscillations: 

These oscillations occurs due to dynamic interaction between channel and the 
compressible volume, which is initiated by flow excursions. 

(c) Acoustic oscillations: These are high frequency (10-100 Hz) oscillations related 
to time required for pressure wave propagation; occurs due to resonance of 
pressure waves. 

(d) Thermal oscillations: These oscillations occur in film boiling due to interaction of 
variable heat transfer coefficients with flow dynamics. 

(e) BWR instability; It is strong only for small fuel time constant and under low 
pressures; occurs due to interaction of void reactivity coupling with the flow 
dynamics and heat transfer. 

(f) Parallel channel instability; there are various mode of flow redistribution; which 
occurs due to interaction between small number of parallel channels. 

However most common among these instabilities is density wave oscillations. 
Present work analyses this instability. If a flow perturbation w' is introduced at the 
inlet to the boiler tube, there will be perturbation in the single-phase pressure drop 
Aj?,'_,i,in the subcooled region, which is in phase with the flow perturbation. There will 
however be a lag in the movement of the boiling boundary, which in turn leads to an 
increasing lag in perturbation in the density and velocity of two-phase region. This 
leads to a lag in two-phase pressure drop perturbation Ap 2 _^ . If the lag in two-phase 

t 

pressure drop perturbation at a particular frequency is 180°, the two- 
perturbations (Api^ and Ap 2 .^) cancel out leaving same constant pressure drop I 
with an oscillatory flow. This is the threshold of instability. When the perturbation ! 
crosses this threshold then, the two-phase pressure drop, at this particular frequency i 



lags by180°, and has amplitude greater than that in single-phase region, the overall 
restriction of constant pressure imposed by the headers will lead to divergent flow 
oscillations. 

Rizwan-Uddin [2] identified that the axial variation in the mixture velocity is the 
origin of density wave instability in two-phase flow. 

It can be now said that the mechanism behind this type of oscillations is 
now well understood in addition to many other type of instabilities, which are cited 
earlier. However experimental and analytical Investigations are still being carried out 
(a) to predict the onset of instability because threshold of instability and stability 
behavior is found to vary by differences in boiling loop geometry and their operating 
conditions, and (b) to understand the characteristics of oscillations under unstable 
operations. 

1.2 Experimental investigations 

1.2.1 In natural circulation loops 

Experimental investigations in two-phase natural circulation loops having 
single heated boiling channel have been carried out by several workers[3-7 J.They 
observed density wave instability, which was found to increase with increase In 
channel exit restriction, inlet subcooling and decrease in system pressure, channel 
inlet restriction and down comer level. 

Bergles et. al. [3] has carried out experimental investigations in a natural circulation 
loop similar to thermosyphon reboiler used for chemical processing and waste heat 
recovery. They also observed that when the liquid level in downcomer is above the 
riser exit, the system is stabilised. Lee and Ishii [4] found that the non-equilibrium 
between the two phases such as during flashing, created flow instability in the loop. 
Kyung and Lee [5] investigated the flow characteristics in an open two-phase natural 
circulation loop using freon as the test fluid. They observed three different modes of 
oscillation with increase in heat flux; 



(a) Periodic flow oscillations characterized by an incubation period of boiling 

(b) Continuous circulation, which is maintained with churn/wispy-annular flow 
pattern. This was found to be a stable operation mode in which the flow increases 
with heat flux and then decreases with further increase in heat flux. 

(c) Periodic flow oscillations characterised by continuous boiling inside the heater 
section, (i.e. there is no incubation period); and void fluctuates between 0.6 to 1.0 
repeatedly. In this mode mean circulation rate was found to decrease with increase 
in heat flux although mean void kept on increasing. 

Jiang et.al. [6] Observed three different kinds of instabilities- 
geysering, flashing and density wave oscillations during the start-up of natural 
circulation loop. Wu.et.al. [7] observed that flow oscillatory behavior was dependent 
on heating power and inlet subcooling. Depending on the operating conditions these 
oscillations can be periodic or chaotic. Mathisen [8] has carried out experimental 
investigation in natural circulation loop having parallel heated boiling channel. He 
observed that increase of input power reduces the flow stability. He also found that a 
possible remedy is an increase of system pressure to stabilize the flow at high input 
power. Fukuda and Kobri [9] have classified density, wave instability as type I and 
type II depending on the operating condition of the system.Type I instability occurs 
under low quality condition due to the dominance of gravitational pressure loss term. 
Type II instability occurs under high quality condition due to the dominance of 
frictional pressure loss term. Aritomi et.al. [10] observed three kinds of instabilities 
during the power raising process of a natural circulation loop with twin boiling 
channels, such as geysering, in-phase natural circulation oscillation and out-of- 
phase density wave oscillations. Saha et. al. [1 1] observed steady state and stability 
behavior of a natural circulation loop. They observed that geysering like instability 
was occurring during start up which developed into type I instability with increase in 
power. They also observed that system stability increases with increase in inlet loss 
coefficients and with decrease in inlet subcooling. 



1.2.2 In boiling forced circulation loops 


Similar experiments were carried out to study the flow instability in forced circulation 
loop. Aritomi et. al. [12] found that the flow oscillations in parallel channels were 
independent of magnitude and nature of disturbance at the inlet . They also found 
that when flow conditions differed between the channels (but where the individual 
characteristic frequencies and inertia masses were approximately equal) the 
instability behavior almost agreed with that of under the average operating 
conditions. The system was found to be more stable with increase in difference of 
flow conditions between the channels. Takitani and Takemura [13] observed that 
the phase difference between inlet and outlet flow rates in the boiling channel was 
180 degree. They also observed that threshold of instability was unaffected by 
presence or absence of superheated steam. Nakanishi et. al. [14] observed both 
fundamental and higher mode of oscillations in their experiment. Wang et. al. [15] 
observed density wave oscillations up to 100 bar in their experiments. Xiao et. al. 
[1 6] observed pressure drop and thermal oscillations in addition to the density wave 
oscillations which occurred both at low and high flow quality in their experiment. 
Their experiment also revealed that density wave oscillations can appear up to 1 92 
bar and disappear above 207 bar. Fukuda et. al. [17] conducted experiments in 
multiple parallel channels (upto 7 channels) with short or long riser pipes. They 
observed type 1 density wave instability in both the cases. Collins and Gacesa [18] 
tested the effect of mass velocity and power on flow oscillation frequency in a 1 9 rod 
bundle with steam-water flow at 800 psia. They found that oscillation frequency 
increases with increase in mass velocity as well as increase in power input to the 
channel. Dijkman[1 9] had tested the effect of cosine heat flux distribution. He found 
that cosine heat flux distribution stabilizes the flow. 

1.3 Theoretical investigations in boiling systems 

Several studies based on linear and non-linear analysis have 
been carried out in the past to investigate the flow instability in the two-phase flow 
systems. These analyses solve the energy and continuity equations to obtain the 



enthalpy distribution in the heated region. The momentum equation is then 
integrated around the loop with this enthalpy distribution. In the linear analytical 
technique the governing equations are linearised by perturbing equations around the 
steady state and solved analytically to obtain the characteristic equation and the 
stability of system is investigated from the root of the characteristic equation. In the 
non-linear analysis, the governing equations are solved numerically by using finite 
difference method .The linear analysis takes less CPU time to evaluate the threshold 
of instability. Hence, it is useful when analysing the mechanism of instability by 
plotting stability maps. But the accuracy of this method is limited. Beyond the 
threshold, non-linear oscillations may appear which can not be predicted by the 
linear stability technique. However, the problem with non-linear analysis is the 
numerical diffusion which may appear in the calculations and then it is difficult to 
interpret the physical instability of the system. 

1.3.1 Linear analysis 

Investigations based on linear stability analysis for natural circulation loops have 
been carried out in [20-25]. They assumed homogeneous two-phase flow in their 
analyses and validated their predictions with the test data. Ishii and Zuber’s[20] 
model was modified by Saha and Zuber [21] by taking into account the thermal 
non-equilibrium effect between the phases. They found that thermal non-equilibrium 
effect between the phases predicts a more stable system at low subcooling when 
compared with thermal equilibrium model. Furutera [22] found that the threshold of 
instability depends on the two-phase flow friction multiplier and the heat capacity in 
the subcooled boiling region .Lee and Lee [23] predicted the threshold of instability 
for Ledinegg and density-wave oscillations and showed that the density wave 
instability analysis is sufficient for the stability analysis of the two-phase flow system. 
Wang et al. [24] predicted the threshold of density wave instability for a two phase 
natural circulation loop and found that instability can appear at low as well as high 
power levels. The stability of density wave oscillations was found to decrease with 
decrease in (a) diameter of riser pipes and (b) system pressure. Nayak et.al. [25] 
performed the linear stability analysis to predict the threshold of instability for a 
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typical configuration of AHWR. They found that density wave threshold decreases 
with an increase in pressure and a decrease in inlet subcooling. They also found that 
existence of power distribution among the channels reduces the density wave 
instability if the channel exit quality is same. 

Similar analysis has also been carried out by numerous investigators for 
forced convection boiling systems [26-29]. Sumida and Kawai [26] derived 
expressions for the hydraulic stability of the boiling channels in terms of flow 
impedance (defined by the dynamic response of pressure drop to the inlet flow). 
Based on these they classified the density wave instability into three main 
categories. Fukuda and Hasegawa [27] derived the mode of oscillations of parallel 
boiling channels from the characterstic equation. Taleyarkhan et. al [28] considered 
subcooled boiling, heat flux distribution in the channel, heater wall dynamics, slip 
flow and transient mixing of flow paths in their analysis. They found that subcooled 
boiling destablises the system and the system stability improves with ventilation 
between the channels. Clausse et. al. [29] found that for system having identical 
channels(a) system can oscillate with all the channels in-phase with the external 
loop (b) channels themselves oscillating out of phase maintain a constant flow rate 
in the loop. If the channels are different, complicated modes of oscillations involving 
coupling between the channels and loop can be formed. 

1.3.2 Non linear analysis 

Gurugenci et. al [30] developed a numerical code to generate limit cycles of 
pressure drop and density wave oscillations in a boiling up-flow system in a channel. 
Chatoorgoon et. al. [31] developed a computer code named as ‘SPORTS’ which 
solves the conservation equations numerically with minimum approximations, 
avoiding the use of property derivatives and matrix inversions. It also permits large 
and small time steps. This code can predict the limit cycle oscillations in both natural 
and forced circulation systems. Nigamutlin et. al [32] found that the heat transfer 
coefficient does not effect the threshold of stability while the thermal wall inertia has 
a strong stabilising effect. 



1 .3.3 Relevant to power reactors 


Apart from all these fundamental studies cited above, other investigations exist on 
flow instabilities in power reactors. Krishnan and Gulshani [33] investigated the flow 
instabilities in the CANDU type PHWRs analytically and experimentally under two- 
phase parallel channel natural circulation conditions. They observed intermittent flow 
oscillations that occurred at low exit qualities caused by the periodic formation and 
expulsion of slug bubbles in the boiling channel: and, an oscillatory flow instability 
that occurs at high exit quality characterised by large amplitude sinusoidal type flow 
oscillations in the channel. Mochizuki [34] studied the instability behavior of the 
heavy water pressure tube type reactor in a simulated experimental facility both 
under natural and forced circulation conditions. Also, he carried out numerical 
investigations to clarify the various types of flow instabilities in the reactors. 

All these studies clarified that the density-wave oscillation is 
the most common type of instability that can occur both in boiling natural circulation 
loops and in forced convection systems. Occurrence of this instability depends on 
the system pressure, channel inlet and outlet resistances, heater power and inlet 
subcooling conditions. 

1 .4 Objective of this thesis 

For two-phase flow instabilities more accurate predictions are being sought. 
Eariier most of the codes developed for this purpose are based on the 
homogeneous model. In homogeneous model relative velocity between vapor and 
liquid has not been considered. However, in two phase flow there is always a 
relative motion between the liquid and vapor phase. In drift-flux model relative 
velocity between liquid and vapor is considered by taking into account the vapor drift 
velocity (Vgj) in the two-phase region. Earlier some codes based on drift flux model 
have been developed. But these codes are applicable to forced circulation loop only. 
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The objective of this thesis is to develop and validate the code 
based on drift flux model for prediction of stability limits and to study the effect of 
various parameters on stability. 
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CHAPTER 2 


DRIFT FLUX MODEL 

2.1 Modelling Of The Stability Behaviour of Two-phase Natural 
Circulation Loop 

The mathematical model to be presented in this section has been applied to 
investigate the stability behaviour of single channel two-phase natural circulation loop such 
as shown in Fig. 2.1. Relative velocity between vapor and liquid is included by using the 
drift-flux model for two-phase region. In the analysis it is assumed that fluid is 
incompressible in the single phase region. 

In the two-phase region we assume that 

(a) the two-phases are in thermodynamic equilibrium, 

ps 

(b) pressure profile in the loop is always in a quasi steady state i.e. — = 0, 

dt 

(c) no carry-over or carry-under in steam drum (i.e. complete separation), 

(d) complete instantaneous mixing of feed-water in the steam drum, 

(e) heat losses in the pipings are negligible, 

(f) properties don’t change significantly along the loop, since Ap across the different 
components (riser, downcomer, cooler, heater and fittings) is « than system 
pressure, 

(g) inlet subcooling ATsub is constant for a particular power and feed-water temperature 
condition. 
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Fig 2. 1 Schematic of natural circulation loop considered in the modelling 




With these assumptions the conservation equations of mass, energy and momentum 
for one-dimensional two-phase flow are given by 


Conservation of mass of mixture 


— +—(PmK) = ^ 

ot oz 


Conservation of mass of vapor phase 




Conservation of energy of the mixture 


^ . r/ ^ ( Pf ~P- PsPt 


o — "^+0 V — ^ = LJLL.11 g y h 

St " az A, az p, t^p 


Conservation of momentum of the mixture 


, = P l^ + v ^1+Ap v^+p g + 

/^m\ ^ 5z J 2D ^ ^ j ^ 


d Pf Pm P fPg 
^\p.-Ps Pm 


These equations have been taken from Belblidia et. al.[38] 

The thermal equation of state gives the density as a function of saturation pressure 

Pf=Pf(Ps); Pg=Pg(ps) pm = apg+(l-a)pf (5) 

Caloric equation of state can be written as: 

hf = hf (ps); hg = hg(ps) hm = {apghg + ( 1 -a) pfhf} /pm (6) 

The vapor drift veIocity(T'^.) is velocity of vapor phase with respect to velocity of center of 


volume or volumetric flux density of the mixture. 
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The volumetric flux density of the mixture is given by 


j = {\-a)Vf+aV^ (A.5) 

so vapour drift velocity is given by 

''zi-K-J 

r„=(i-aXn) (A.9) 


Where is relative velocity of vapour with respect to liquid phase. 

So it accounts for the relative motion between two phases. If Vg,- = 0.0, then present model 
will be equivalent to Homogeneous model. 

2.2 Steady State Equation 


The governing equations for the steady state conditions are obtained by dropping the 
time derivatives Ifrom equations (1) to (4). 

2.2.1 For single-phase region 


Conservation of mass equation becomes 


dz 


{PfVf) = Q 


Conservation of energy equation becomes 

PfV r — for heated region 
dz 

=0 for adiabatic region 

Conservation of momentum equation for the liquid phase 


(V) 

(a) 

(b) (8) 



(9) 


-Sp 

dz 


-PfVf 


dVf f 2 
dz ID ^ ^ 


^PfZ 


acceleration viscousloss gravitaional 
term term term 


2.2.2 For two-phase region 

Conservation of mass of two-phase mixture 

^(.p.K) = o 

dz 


Conservation of vapour mass for two-phase becomes 

/ \ 




^PfPz y 
V Pn, 


for heated region 


dz 


( 

— 

V Pm 


^PfPz y 


for adiabatic region 


Which can be transformed to (refer equation B. 1 1 to B.29) 

for heated region 
dz 

— 0 for adiabatic region 

Conservation of energy of the mixture 


p V 

r'm n 


_ qlPh ^ 

dz dz 


Pf-PmPfP, y. 

I A;, 


for heated region 


dz 


^ Pf- Pm PfP, y. \ 
\ Pm A/7 ® 


for adiabatic region 


Ck=j + Vgj (This is an approximation for slug flow) 

Pm n ^ 

Pm 


Q. 


where Fj, =vapour generation rate= 


PfP> 


( 10 ) 

(a) 

(b) (II) 

( 12 ) 

(a) 

(b) (13) 

(A.14) 
(A. 10) 
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Mf, PfPz 

Conservation of two-phase momentum 


— dp dV„ f„ ^-2 d 

=■ O 1 P Q PL 

oz oz zD oz 


P f Pm PfPg 
Pm- P i, Pm ^ 


\ 


( 14 ) 


J 


acceleration viscousloss gravitaional slip loss term 
term term term 


2.3 Steady state solution 

2.3 (a) For the single phase region 

Equations (7) to (9) for the single-phase region can be solved analytically and the expressions 
for pressure drop for single-phase region can be obtained as 

K,PfV,J 


(Ap)i-^ = pfS^ + 


fsPfVfisiL 

2 D 


+ - 




ID 


- for adiabatic region (a) 
for heated region (b) (15) 


In equation (15) acceleration term is zero due to conservation of mass and incompressibility 
of fluid. 

A-ss = length of single phase region in heated channel 
AZ = change in height (elevation difference) 


2.3 (b) For the two phase region 


Equations 10-14 can be solved analytically to obtain the pressure drop for 2-(j) region. 

aFL 


£ dV 


^l^PmsiKss^^mss (SinCC P P f,J f,,, ) 
= o V f dV =n V \v -V 1 

r mss mss ^ss r mss mss L mss^L mss,^s J 

(2) viscous loss term= 


2 D 
fm 

p„. 

Now since p V ^OrVr 

r mss^ mss H fm fiss 


f 

JAs: 


dz 
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Viscous lossterni= — - — —dz t ^ 

J<“2B p„, 2D^ ' 

XT ,. P fiss^kss,X^ 

Now from equation 

^kss 

=^p^f,s.Vl, f ^ dz 

now 6z= -^ — ^ 


(B.40) 


so VISCOUS 


( from equation B.36) 

loss term f = /-l 

on r ks a,, AD ^ 


kss,Xss 


gravitational loss term= PmssSdz 


using /?„„ =ED^JsLh^ ^we get 

^kss,L 


^hss 


dz 


• . dC^, 

using dz = — — 


(B.36) 


so gravitational loss term = ^ [in 


InC 


kss,Ass 


(4) Slip loss term= £ 


\ 


dz 


Pf Pm PfPs j/2 
Pm P g Pm j 


P f Pm P fPg |-/2 
\Pm-pg Pm 




P [ Pmss,L P fPg yl 

^ si 


Ms ^ 


PmssX-Pg Pd 


J 
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/A \ T \t 7 T/ 1_l P Jiss^ fiss l»l» I Tfi^tn 

(Ap)2-(D = Pm,s^mssWmss,L 7^~ ^ ^ kssMs ^ 


4D Cl C 

H-zy ^^sS^kss,Xss 


P fssS^kss.X^ 


4- 

Pf PmssX 


* Cl 

ss 


^ P mss P g 

Pm 


for heated region (a) 


7 

and = P„.sS^+^Pn.sVLL^ ^'—^-- for adiabatic region (b) (16) 


It is important to note that AZ is elevation difference in the direction of flow. If the flow is 
from higher to lower elevation then AZ will be with a negative sign. 

now (Ap)total ~ 2j(Ap)i.(j)adiabatic “*'2j(Ap)i-<j)heated 2!!(Ap)2Hj)heated ‘^2(Ap)2-(()adiabatic (f 

since it is close loop with no pump. 

(^p)total=0 0 

So I(Ap), -<|)adiabatic "^2(Ap)i-({)heated 2(Ap)2-<jiheated "^2!!(Ap)2-ij)adiabatic 0-0 (1^) 


For a assumed flow rate value of pressure drop can be calculated, since all other variables 
such as Ck, Vm, Vf, can be obtained. So equation (18) can be solved iteratively to get the 
steady state flow rate. 

2.4 Linear Stability Analysis 


2.4.1 Perturbation Variables 


The set of conservation equations (1) to (4) are linearised by superimposing of small 
perturbations of , W' , , h'j., p'„ , C[ over the steady state such as follows: 


where f; = , F^ = V^se^‘ , h} = h^ee^‘ p'^ = p„£e^\ Ck' = Ckse^' (20) 

Vj;,hf, V^,p^, Q are average of perturbed variables and s is the stability parameter. 


The governing (conservation) equations are linearised and solved analytically. Kinematics of 
the flow is solved first. Then the momentum equation is integrated to obtain the response in 
total pressure drop perturbation to an inlet flow rate perturbation. 

2.4.2 Kinematics of the flow 


2.4.2 (a) For the single phase heated region 

The perturbed energy conservation equation for single phase heated region can be 
written as (refer equation B.l to B.5) 


dh) 

dz 


+ 


sh', -q:p„v; 

^fiss Pj^flSS^ 


( 21 ) 


Solution of equation (21) is 
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h' 

PfAs A‘ss 
from equation (8) 


^SS 
V 


1-e > (refer equation B.6 to B.IO) 

CP,, 

pfVf,ssA 


( 22 ) 


(23) 


For a uniform axial heat flux, the enthalpy increases linearly with Z. A positive perturbation 
in enthalpy at the boiling boundary causes a negative perturbation in the boiling boundary. 
Therefore, this can be expressed as 


dh 




^fss,Xss 


X' = - 


PfAVj 




q'lA 


fiss 


h! 


^fyZSS 


or 


_ P f h ? w^h 

Pf^c^ ^fiss 


sZss 


1 — 

hence X' = 

s 

where xi = A-ssA^ fiss is residence time of fluid in single phase region 

2.4.2(b)For 2-0 region in the core 

Volumetric flux equation (refer equation B.l 1 to equation B.20) 

dz 

Density propagation equation(refer equationB.21 to equation B.28) 


(24) 


Spn, I Q Spm 


-Pm^ 


dt dz 

solution of volumetric flux equation is 

j = VfXt)+a{z-i(t)} 

where X{t) is height of boiling boundary from the inlet of heated channel. 
The perturbed portion of kinematic velocity can be written as 
Since Ck ~Vfi+Vgj 

From equation (25) Ci<= F^(t) + Q.{Z -X(i)} + 

On perturbing both sides, we get 


(25) 
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( 26 ) 


So c;=f;-qz' 

c; = F;{i-Qr,(^)} 


c; =T,(sw; 

The perturbed density propagation equation for core can be written as 
(Refer equation B.28 to equation B.34) 

n' — ^k^ssPmss 

dz Ckl 


This can be integrated from to Z using the fact that a positive perturbation in the boiling 
boundary causes a negative perturbation in mixture density. To a first approximation it can 
be written as 


dz 


mss ) 


Q, 


Pf 


C 

^kss,^j;s 



(27) 


o' - 1' 

Pm, Ass - 


Ck. 


'ss,Xss 


Thus, the solution for the mixture density is as follows(refer equation B.35 to B.48): 




Pm =Pf 




V 


kss,Ass 


\^~^ssK ^kss J 


c 


1 


■C 


"'ksSyXss 




n !'■ 


Q. 


s-Q. 


c: 


ss J 


pf 


V 




\^ksSyAss ^ksSyL J 


c. 


ksSyXss 


(28) 


On perturbing equation (B.17) and linearising, the solutions for the mixture velocity 
perturbations are obtained.(refer equation B.49 to B.51) 


V’„=T,{s)V'+- 


In 

Pf 


Ck, 


V, ^^ss,Ass J 


(29) 


2.4.2(c) For unheated riser section 

Since the riser is unheated, Q, = 0.0 
Volumetric flux equation 

dj 


dz 


^0 


On perturbing above equation and solving 


( 30 ) 
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Ck = j' = constant 
density propagation equation 

+ Q = 0 

dt dz 

On perturbing equation (B.52) and linearising(refer equation B.52 to B.55) 

^ Pm + ^kss ~ ~ ^ 

oz 

On solving equation (37) density perturbation is obtained: 


Pm ~ Pm,l^ 




( 31 ) 

(B.52) 


(B.36) 


(B.55) 


On perturbing equation (B.17) and solving for velocity perturbation 
.(refer equation B. 58 to B. 61) 


V:=C[+p,V^p'„^,e 




(3.61) 


2.4.3Dynamics of the Flow 

By perturbing equation (4) perturbed momentum equation is obtained. Now this equations is 
integrated to obtain the perturbed pressure drop . 


ApL/ ^ Ap; ^ Ap; ^ Apl ^ Ap'^ ^ Ap; ^ Ap'f 
K-o Wj'-cD K-® 


(32) 


where = Pj-A (33) 

2.4.3(a)For single phase adiabatic region 

(i) Inertia component (refer equation B. 79 to B. 80) 

Ap,- _ 

(ii) Acceleration component (refer equation B.65 to B.68) 


(iii) Frictional component (refer equation B.69 to B.72) 

^P'f fsVfissL 

DA 


22 



(iv) Loss component (refer equation B.73 to B.76) 

AP/' _ 


A 


(v) Gravitational term (refer equation B.62 to B.64) 

0.0 






(vi) Slip component (refer equation B.77 to B.78) 

Ap; 


= 0.0 


w, 


2.4.3 (b) For single-phase heated region 

(i)Inertia component (refer equation B.95) 
Ap' _ sL 


wL A, 


(ii) Acceleration component (refer equation B.89 to B.90) 

Ap; 


= 0.0 


W 




(iii) Frictional component (refer equation B.86 to B.88) 




2DA. 


(iv) Loss component (refer equation B.91 to B.92) 

A/?; _ 


w, 




(v) Gravitational component (refer equation B.81 to B.85) 

Ap; _ gr,(.y) 




(vi) Slip component (refer equation B.93 to B.94) 

Ap; 


0.0 


w 


i-d) 
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2.4.3 (c) For 2-phase heated region 

For two-phase region following equations are used frequently 


P'm,L [ 


1 r(s) 


^ ) I ^^ss.L J ^Ks,Xss 


\ Q 1 

£2„r,W - 

A^\ s-Li„ nc. 


1 Y ^kss.Ass 1 ^ 


kss^Xss Jl ^kss^L 


K, r,(s) , vj ^ />:, 


<o pf^ 

(i)Inertia component (refer equation B.122 to B.126) 


In 

J ^ss^c Ois.Asj K -^c ^ ~ ^ss ^kssj 


^ ^ss ^kss,Xss 


+ -f- r,(^)Q„ 
I 


r,wn„'i 1 ( c, 


^is ^ I 


(ii)Acceleration component (refer equation B.127 to B.143) 

Ap: . (, ^ r,(j) f, 


4 I S-^ss) Ws 


^SS ¥(^ fP 

^ 1 V'-^ife5s,I/ \^kss,Xss/ 


f^s -S^ + ^SJ ^ z-’ T ^kss.L 

7~ 775 ^ O 2 W ^ 

\ ^ /L*^ “ ^kss^^ss 

s 

Q ^ c"“ fr \asi: 

^fiss ^ (*^0 ^kss , L ^ ^fiss ^gj ^ 

*^r ^kss,Axs ^kss.Ass ^1-0 


'^hs^Xss / ^ 
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(iii)Frictional component (refer equation B.lOl to B.l 17) 


,^Vl± = Jji_ y r ^ 

DA, y,DA, Cl^ 




V 

2 fiss 




\~i 




+ ■ 


/. 


ID pfA^ 
P Pf 


1-^ 


\^r p\ 


^kss,Xss J 




/7V^r,(5) 

2£>A 


Where . — 


y, 


+ 




ss 

\^^ss ~^J 


a. 


,'Fhizi 
' Cl. 


^ss J 


^ / V±h: 


20„-j 


(iv) Local loss component (refer equation B.l 44 to B.l 47) 




r 


w, 


\-6 


W, 




jy ^ T/2 

, ^mss.L 

V 






(v) Gravitational component (refer equation B.96 to B.l 00) 


Ap; 


8 _ 


a 


w 




sA^ 


r,(s)-r,{s) 


1 




+ • 


1 Q 


ksSyXss 


s — Cl„ A,. 


r,(^) 


kss.Ass 

1 ^ 


^r,(^) 


g 


^kss.^s ^kss,L j 


g- 


(vi) Slip component (refer equation B.l 18 to B.l 21 ) 


^P\ P'm,L Pf Pf^a 




2/?y 

P mss^L ) 
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2.4.3 (d) For riser section 

Following equations are frequently used. Value is given by 

Pm,L, _ Pm,L ^-s{L,-L)/ 

K-O 


( 36 ) 


Kj, _r,(5) 


W,. 








‘'l-® Pf'^c 

(i) Inertia component(refer equation B.161 to B.164) 


Ap,' 


w,_ 


pfA 


m-L)+ 


\ P mss, L J 


r' m,L 






( Qjj,i 1 


exp 


(37) 


(ii) Acceleration component (refer equation A1.165 to B.172) 


Ap: 


w,_ 




s(Lx-l) 


-1 


EEE—Q y 

^ kss,L^ mss,L^ 






mss,L 




(iii) Frictional component (refer equation B.152 to 155) 




w, 




D 


Pmss^mss 


rAs)m -L) 

pfA 


■PA 






w, 


l-<zi 


J 


i,z£ (A - - . a -i.o 

c 

^kss,l 


+ 


2D ”” 




P. 


m,L 


V J 


Q 


kss 


s J 






kss 


(iv) Loss component (refer equation B. 173 to B.174) 
— = 2K,p^, F ^ + K, ^ ri 




\-6 




W, 


1-^ 


(v) Gravitational component(refer equation B. 148 to B. 151) 


Ap; 






W, 


A- 


r-c. 






V 


ksSjL 

\ ^ y 


;exp^^l— ^-1.01 

^kss.L 


(vi) Slip component (refer equation B.156 to B.160) 


ap: -'2-vip}p, 


w, 




Pmss 


P m,Li 

\ 


P 


m,L 


W,. 


+ 


PfPsK 


2 / 


2 

Pmss 


P m,ly Pm,L 






J 


Now for the loop can be calculated. 
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N?,o,ai = {Ap„ L, + (Ap'L«/ = 0.0 for a closed natural circulation loop 
but {l^p^Xotai ~ 0-0 from equation (18) 

so(Ap'L/ =0-0 

but =E(Ap')i^,adiabatic +Z(Ap')l r<j)heated E(Ap )2Hjiheated "^ECAp )2Hjiadiabatic 0.0 (38) 

Equation (38) is the characteristic equation. It is solved for the stability parameter (complex 
no.).If the s is the root of the characteristic equation given by 

S = 7] + jtJJ 

Where rj is the real part and cr is the imaginary part of the root. Then the system is 
considered to be stable if all 7 < 0 , and unstable if any 7 > 0 .at the neutral threshold point at 
least one of the 77 =0. 

Decay ratio defined as the ratio of the successive amplitude of impulses can be 
obtained as: 

2^7J 

DR=exp''^^ (39) 

The decay ratio indicates the stability margin, which the system can have at, any operating 
condition. If decay ratio is less than one then the system is stable. If it is more than one the 
system is considered unstable. If it is equal to one then the system is at threshold of stability. 
In this condition system is called neutrally stable. 
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CHAPTER 3 


MODEL VALIDATION 

3.1 Experimental Loops 

It is required to validate the theoretical model discussed earlier in chapter 2. For this 
purpose the code TINFLOS-D has been developed and validated with experimental data 
obtained from two natural circulation loops in BARC- (a) HPNCL (b) Apsara.The first 
experimental loop HPNCL is in Hall No. 7 of BARC.The second natural circulation loop 
Apsara is looted near Apsara reactor in BARC. These loops along with the geometry are 
shown in FigS.l and Fig 3.2 respectively. 

Code TINFLOS-D has also been validated with two other natural circulation loops reported 
in literature. One experimental result obtained from Chexal and Bergles [3] and other 
experimental data obtained from M. Furutera [21]. These loops and their geometry are 
shown on Fig 3.3 and Fig 3.4 respectively 

3.2 Validation 

Fig. 3.5 shows the comparison of predicted mass flow rate with experimental mass flow 
rate for HPNCL. It can be seen that prediction are in good agreement with test data. Fig 3.6 
and Fig 3.7 also shows the threshold power of stability for experimental data and 
TINFLOS-D at 1°C and 2.5°C subcooling respectively for the HPNCL. Table 3.7 also 
shows that predictions are in agreement with the experimental results for HPNCL. 

Table3.8 shows the validation of predicted data with experimental data 
for Apsara loop. It can be seen that predictions are in good agreement with the experimental 
results. Fig. 3.8 and Fig.3.10 show validation of predicted data with experimental data for 
Furutera’ s loop and Bergles loop respectively. It can be clearly seen from these figures that 
predicted data are in good agreement with the test data. 
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In Furutera’s loop there are two Spherical buffer volumes as shown in Fig 3. 3. these 
spherical volumes will act as dampener for oscillations. But in present analysis the effect of 
these buffer volumes has not been considered, and these buffer volumes has been 
considered as the pipe of constant diameter 0.0215 meter. Diameter of Steam separator is 
not available in literature. So it has been considered to be of very large diameter, so that 


area ratio 


A 

A 


« 0. Value of throttling loss coefficient and the orifice loss coefficient is 


not available. So hit and trial procedure has been adopted for proposed throttling and orifice 
loss coefficient. Fig. 3.9 shows the comparison of experimental results with the predictions 
by 7 models .Fig. 3.9 and table 3.6 are reported in Furutera et. al.[21].It is given in the 
literature that the prediction by model DST-1 is best followed by DST-0 and DST-6.on 
comparing the errors (table ) we find that Prediction by DST-1 is most accurate followed by 
TINFLOS-D(Vg,_su.g), DST-6, TrNFLOS-D(Vg,-0.0), DST-0, DST4,DST-2,DST5 and DST- 
3 respectively.DST-1 is predicting better than the TINFLOS-D(Vgj_siug).It is due to the fact 
that in DST-1 subcooled boiling is considered, while in TINFLOS-D(Vg]_s]ug) subcooled 
boiling has been neglected. 


In Bergles loop there is a scope for improvement in simulation, since subcooled 
boiling has not been considered. Also in the Bergles loop value of inlet valve loss 
coefficient and orifice loss coefficient is not available in the literature. So hit and trial 
procedure has been adopted for the proposed inlet valve and orifice loss coefficient. 


Table 3.1 Errors associated with predictions 


Model 

sum of error square/N 

DST-0 

291.711 

DST-1 

11.86 

DST-6 

238.425 

TINFLOS-D(Vgj_s,ug) 

44.134 

TINFLOS-D(Vg,=0) 

259.87 
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Region 3 

All dimensions in meter 

Fig 3.2 Apsara loop 
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•^2 


2.077 




AM dimensions in inch 

3.4 Bergles loop 



Table 3.2: Data for HPNCL Loop 



Regions 


1 

2 

3 

4 ' 

5 

6 { core ) 

7 

8 

9 

Az 


loom 




■nH 




L 






1.18 

0.405 

1.11 

fBSM 

Klexit 

0.4526 

1.185 

1.185 



0.0 

0.102 

1.185 

0.767 


Table 3.3: Data for Apsara Loop 



Regions 


1 

2 

3 

4 

5 ( core ) 

6 

7 

Az 

- 0.25 

- 2.227 

0.0 

0.39 

0.29 

1.797 

0.0 

L 

0.25 

2.227 

2.10 

0.39 

0.29 

1.797 

2.07 

Klexit 

0.4936 

1.185 

1.185 

0.0 

0.0 

1.185 

0.9663 


Table 3.4: Data for Furutera’s Loop 



Regions 

1 

2 

3 

4 

5 

6 

7 

lEBlSi 

9 

10 

Az 

- 0.30 

- 1.7655 

0.0 

0.0 

- 0.3115 

0.0 

0.3115 

1.831 

0.2359 

0.0 

L 

0.30 

1.7655 

0.01 

0.034 

0.3115 










1.185 

1.185 

1.185 

0.0 

0.0 

1 . 185 

1.0 


Table 3.5: Data forBergles Loop 



Regions 


1 

2 

3 


5 

6 

Az 

- 1.0795 


0.15875 


0.15875 

0.0 

L 

1.0795 



moogjii 


0.5334 

Klexit 

0.495 

1.185 

0.0 

0.0 

1.185 

0.97 
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Fig 3.5 Comparison of predicted mass flow rate with 
experimental mass flow rate for HPNCL under steady state 
condition 
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012345678 

Pressure(bar) 


Fig 3.6 Threshold power of stability-prediction and experimental data for 
HPNCL 
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Pressure(bar) 
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Fig 3.7 Threshold power of stability- prediction and experimental data 
for HPNCL 
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10 20 30 

Pressure(bar} 


Fig 3.8 Threshold power of stability-prediction and experimental data 
for Furutera’s loop 


38 




Pressure MPa 


Fig. S'.? Comparison between experimental results and theory 


Table 3 '<$ 

Contents of the models 


Mode|\..^ 

Heat capacity 

Fricticnal coefficient 

Curve 

Single -phase 
liquid region 

Subcooled 
boilirq region 

DSTO 

ignored 

ignored 

constant 


DST 1 

considered 

considered 

constant 


DST2 

ignored 

ignored 

transient 

— 

DST3 

considered 

considered 

transient 

— 

DST4 

considered 

ignored 

transient 


DST5 

considered 

considered 
(Jens- Lottes) 

transient 

— 

0ST6 

considered 

ignored 

constant 
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SubcoolingfC) 
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Fig 3.10 Threshold power of stability-prediction and experimental 
data for Chexal and Bergles loop 
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Table 3.7 Comparison of TINFLOS-D prediction with experimental 
observation for HPNCL 



Experimental Conditions / Observations 

Prediction 

Remarks 

Pressure 

Power 

(kW) 

Subcooling 

Stable/ 

Unstable 

Stable/ 

Unstable 

1 

1 atm. 

20.1 

3.55 °C 

Unstable 

Unstable 

Tallies with the 
experiment 

2 

2.67 bar 

22.2 

2.73 °C 

Unstable 

Unstable 

Tallies with the 
experiment 

3 

1.65 bar 

20.0 

3.22 °C 

Unstable 

Unstable 

Tallies with the 
experiment 

4 

5.0 bar 

20.0 

2.7 °C 

Unstable 

Unstable 

Tallies with the 
experiment 

5 

5.33 bar 

24.0 

2.57 °C 

N.S. 

Unstable 

Stable above 27.5 
KW 

6 

8.06 bar 

19.5 

3.2 °C 

Unstable 

Unstable 

Tallies with the 
experiment 

7 

1.18bar 

13.8 

2.69 °C 

Unstable 

Unstable 

Tallies with the 
experiment 

8 

1.39 bar 

20.0 

3.57 °C 

Unstable 

Unstable 

Tallies with the 
experiment 

9 

5.86 bar 

25.2 

2.47 °C 

Stable 

Stable. 

Tallies with the 
experiment 

10 

3.03 bar 

31.6 

2.35 °C 

Stable 

Unstable 

Stable above 34 

KW 


N.S.— NEUTRALLY STABLE POINT 
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Table 3.8 Comparison of TINFLOS-D predictions with experimental 
observation for Apsara loop 


SI. No. 

Power 

(kW) 

SubCoolingCC) 



Prediction 

1 

1.2644 

7.60 

1.504 

Unstable 

Unstable 

2 

1.33736 

4.00 

1.420 

Unstable 

Unstable 

3 

1.3547 

4.00 

1.700 

Unstable 

Unstable 

4 

1.4917 

5.00 

1.980 

Unstable 

Unstable 

5 

1.5886 

8.12 

3.300 

Unstable 

Unstable 

6 

1.6225 

7.32 

4.40 

Unstable 

Unstable 

7 

1.7356 

7.00 

4.760 

Unstable 

Unstable 

8 

1.7576 

4.28 

3.700 

Unstable 

Unstable 

9 

1.8695 

1.90 

4.500 

Unstable 

Unstable 
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CHAPTER 4 


SENSITIVITY ANALYSIS OF HPNCL 

4.1 Prediction of the stability behavior of the loop 

The effect of various geometric and operating conditions on the flow stability behavior of 
the High Pressure Natural Circulation Loop (HPNCL) have been investigated by the 
computer code TINFLOS-D.The geometry of the HPNCL has been shown in Fig 3.1. In 
this HPNCL loop thermodynamic quality of the mixture is very low and is less then 2 
percent. It is interesting to analyse the stability behavior of such a low quality natural 
circulation loop. The sensitivity analysis of HPNCL has been done for lower threshold of 
instability. Upper threshold occurs at higher quality. Since present analysis has been done 
for low quality HPNCL, our results are true for lower threshold only. 

From Fig. 3.5 it can be seen that natural circulation flow rate 
depends on the system pressure. At low pressure, there is increase in flow rate with 
increase in system pressure. The initial increase in flow rate with increase in pressure is 
due to reduction of two-phase friction pressure drop because the two-phase friction 
multiplier reduces with increase in system pressure. At higher pressure, the void fraction 
reduces with increase in pressure and hence the two-phase mixture density increases with 
increase in pressure. This results in reduction of driving head at higher pressure causing 
the flow rate to reduce. 

Density wave instability occurs normally due to the dominance of two-phase 
pressure drop over single-phase drop of the system. This is due to the fact that two phase 
pressure drop is out of phase with the change of inlet flow, while single phase pressure 
drop is in phase with the change of inlet flow. From the figures it can be observed that 
with the increase in power flow stablises. This can be explained from the fact that in low 
quality conditions gravitational pressure drop in two phase is dominant since frictional 
loss is given by: 
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^Pf 


V 




So it is clear that it is inversely proportional to the density of the mixture, it will have 
small value and will not effect the stability characteristic. With increase in power since 
the quality increases, hence the void fraction also increases, since quality is given by the 
expression: 


X = ■ 


1 


1 + 


1 




/ _ A 


--1 




K^sA'sJ 


F 


The gravitational loss reduces due to decrease in two-phase mixture density. Hence the 
flow stablises with increase in power. It can also be observed from the figures that with 
increase in pressure threshold power reduces. This can be explained from the fact that 
with increase in pressure flow rate reduces . Outlet enthalpy is given by 


^out,ss ^inss 




^ P f P mss, out P f P g y , 

Pmss,out A/7 


Where and are inlet and outlet enthalpy to the core. Due to decrease in flow 
rale enthalpy increases. Since the quality is given by; 



So quality also increases. Hence gravitational loss reduces due to decrease in two-phase 
mixture density, which accounts for the reduction in threshold power. 
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4.2 Effect of riser height 


From Fig 4. 1 it can be observed that with the increase of riser height flow stability 
decreases. It can be explained as follows. 

With the increase of riser height downcomer height also increases. So it 
increases the driving force, which results in increased flow rate. Due to increase in flow 
rate quality decreases, as it is clear from the equation mentioned in section 4.1 for quality. 
This increases the gravitational pressure drop. Also the increase in riser height further 
increases the gravitational pressure drop. 


All these effects tend to destabilize the system and the threshold power 
for stable flow increases. 



Fig 4.1 Effect of riser height on threshold of instability 




4.3Effect of downcomer diameter 


The effect of downcomer diameter on the flow instability is 
shown in Figure 4.2 .It is clear that flow destabilizes with increase in downcomer 
diameter. It can be explained as follows 

With the increase in downcomer diameter the 


frictional loss reduces 


APr.. = 


2Dp, 


A‘ 

mss^ J 


since frictional pressure drop is inversely 


proportional to the diameter. It results in increased flow rate. Due to increased flow rate 
two phase flow quality decreases (as it is clear from the equation mentioned in section 4.1 
for quality). It increases the gravitational pressure drop. Hence it destabilizes the flow. 



Fig 4.2 Effect of downcomer diameter on threshold of instability 
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4.4 Effect of core inlet loss coefficient 


It is clear from the Fig 4.4 that flow stabilizes with increase in inlet loss 
coefficient. Due to increase in inlet restriction the pressure drop increases 


^Plss = 


2A^Pf 


. So the flow rate decreases. Hence the quality increases. So the 


gravitational loss reduces due to decrease in two-phase mixture density. So it stabilizes 
the flow. 



Fig4.3 Effect of core inlet coefficient on threshold of instability 


4.5 Effect of core outlet loss coefficient 


It is clear from fig 4.4 that outlet restriction also stabilizes the flow. It can be seen 
from the Figure 4.5 that it is stabilizing little more then the inlet restriction. 


It is due to the fact that in the two-phase region fiictional pressure drop is more due 

( K w 

to' the less density relative to single-phase flow. Since here ' ‘ 


APfo = 




and 


Pxms ^ Pf So here flow is relatively more reduced .it can also be seen from Fig.4.6. So 

quality in this case is higher than that with inlet restriction (Fig 4.7). Hence outlet 
restriction is stabilizing more than inlet restriction. In general outlet orificing is found to 
stabilize less than the outlet orificing. 
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Fig4.4 Effect of core outlet coefficient on threshold of instability 



Fig 4.5 Comparison of effect of core inlet and outlet coefficient on 
threshold of instability 


48 





Flow rate(kg/s) 


1.6 

1.4 

1.2 

1.0 

0.8 

0.6 

0.4 

0.2 

0.0 

0 10 20 30 40 50 60 70 80 

Power(kW) 

Fig 4.6 Comparison of effect of core inlet and outlet coefficient on 
mass flow rate 
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Fig 4.7 Comparison of effect of core inlet and outlet coefficient on 
quality 
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4.6 Effect of inlet subcooling 


It can be seen from Figure 4.8 that with increase in inlet subcooling flow stability 
decreases. Due to subcooling the inlet enthalpy „of the fluid entering the core reduces. 

Hence the outlet enthalpy also decreases. So the quality is lesser (as it is clear from the 
equation given in section 4.1 for the quality). Hence gravitational loss increases. So, it 
tends to destabilize the flow. 



Fig4.8 Effect of intlet subcooling on threshold of instability 


4.7 Effect of vapour drift velocity (Vgj) 

It can be seen from fig 4.9 that with increase in vapour drift velocity flow 
stabilizes. This can be explained by referring Fig 4.10.It is clear from Fig 4.10 that for the 
same change in quality there is less change in void fraction in the flow having higher Vg. 
So for same change in quality there is less chmge in density (p'„) for the flow having 
higher Vg, So the perturbed gravitational pressure drop reduces with increase in Vg. So 
with increase in vapour drift velocity flow stabilizes. 
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Fig4.9 Effect of vapour drift velocity on threshold of instability 



Fig 4.10 comparison of rate of change of void fraction vrith respect to 
quality 


VTT 
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CHAPTER 5 
CONCLUSIONS 

The already existing computer code TESfFLOS has been modified by using the drift flux 
model and computer code TINFLOS-D has been obtained. Model predictions were 
validated using the test data obtained firom the four different loops (HPNCL, Apsara, 
Furutera’s, & Bergles). Parametric sensitivity analysis has been carried out to study the 
flow instability behavior in HPNCL under low quality conditions (i.e. up to 2% quality). 

The following insights have been obtained from this study for lower threshold of 
instability: 

( 1 ) The natural circulation flow rate depends on the system pressure, which initially 

increases then decreases with increase in pressure. 

(2) With increase in system pressure system stability increases. 

(3) The flow instability reduces with increase in power. 

(4) Instability increases with increasing the riser height as evidenced by Fig. 4.2 

(5) Instability increases with increasing the downcomer diameter as evidenced by 
Fig. 4.3. 

(6) Increasing the resistance at the core inlet and outlet increases the stability as 
evidenced by Fig. 4.4 and 4.5. 

(7) The effect of increase in inlet subcooling is to destablise the flow as evidenced 
by Fig.4.9. 

(8) Flow stablises with increasing the vapour drift velocity as evidenced by Fig. 
4.10. 

So, it is recommended that system pressure and power should be high to enhance the 
stability. 

5.1 Suggestions for future work 

In the present code TINFLOS-D the vapour drift velocity has been assumed constant. It is 
taken as the drift velocity for the slug flow. While the vapour drift velocity changes with 
the change of flow patterns. So the vapour drift velocity based on flow pattern should be 
used in future to improve the code TINFLOS-D. 
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The vapour drift velocity for different regions is given below. 

(a) Bubbly regime; [SS] 

(G< 2000 kg/m^ -jand 0.01 <a< 0.01) 

V 1-41 \ scr(p,-p^) Y 

(b) Slug regime: [36] 

(G < 2000 kg/ m ' - 5 and 0.2 < « < 0.65) 

V 0-345 f gA, (a -Pj p 

(c) Annular regime: [3 7] 

(G < 2000 kg/ w ' - 5 and 0.85 < a < 0.90) 

V 

J/ = _2 j 

f /?^(76-75a)y ^ ap^ 

^ p,-sfa j 

(d) Chum-turbulent regime: [37] 

(G> 3000 kg/m^ -s and for all void fractions) 

V 

y - _2 

(C„-l) P. 
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With Co (Void distribution parameter) —1.1 and cc restricted to a maYimiim value of 0.8. 
For values of a > 0.8, {a = 0.8) is used. 

(e) Fora < 0.005 , = 0.0 

(f) Fora = 1 , = 0.0 

(g) Relative velocities in the transitions regions (i.e., 0.005<a <0.01, 0.1 < a <0.2, 
0.65<a<0.85, 0.9<a<l, for G< 2000 kg/m^ -s, and for 2000<G<3000 kg/m^ - 5 ') 
are linearly interpolated in void fraction and/or mass velocity, G. 



Void fraction (%) 

Fig 5. 1 Flow regime map for slip flow correlations 


Sensitivity analysis has been done for lower threshold of instability, which occurs 
at low quality. On increasing the power at a particular pressure, upper threshold of 
instability may be obtained at high qualities (around 15% quality). So in future HPNCL 
should be investigated for upper threshold of instability also. 
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Appendix A 

Identities and definitions with related derivations 

a = local volumetric concentration of mixture (Void fraction) 

a=ElZPj^ (A.1) 

P f ~ Pg 

from p„, = ap^ + (1 - a)p^ (A.2) 

Let Vf, Vg = local velocities of liquid and vapor respectively 
The volumetric flux density of the liquid is 


jf=(l-a)Vf (A-3) 

and the volumetric flux density of vapor is 

j, = aV, (A.4) 

The volumetric flux density of the mixture is 

j = jf + jg = (1-a) Vf + aVg (A-5) 

Defining the mass velocity of the centre of mass of mixture as 

ap,V,+<,l-a)p,V, 


From equation (A.5) 

j=aVg + (l-a)Vf 
= Vf+a(Vg- Vf) 

= [pmVf + pmOt(Vg - Vf)] / pm 

ap/f + (1 - a)pfVf + cc{(XP<. + 0 - cc)p^Wg " 

j _ 
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Pm 


+ a(l - a)PfV^ - a(l - a)p^V^ + ap^V^ 
Pm 


m * ■ 

Pm 


j.y ,<Pf-P,)^~CC)V, 

Pm 

(A.7) 

Where = Vg - Vf 

(A.8) 

V^ = (l-a)W 

(A.9) 

j = V 

J m m 

Pm 

(A.10) 

V = / - ^P y 

p. “ 

(A. 11) 


The kinematic wave velocity is defined as 

dV„. 


C,-J + V^+a 


a 


da 


(A.12) 


However, if Vm does not depend on a. 


dV , 

Then — ^ = 0.0 
da 


In slug flow 0.345 


gP>>,{Pf-Ps) 


Pf 


(see reference [35] ) 


(A.13) 


It is clear from equation (A.13 ) that Vg does not depend on a for slug flow 


So — ^ = 0.0 

da 
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so Ck =j + Vg, 


(A.14) 
(A. 15) 


j Ck - V gj 
Ck=j+Vg-j 

or Ck = Vg =local gas velocity 

Substituting equation (A. 15) into equation (A.1 1), we get 


F =C - V ■ 




P, 


m J 


C-F., 


ap^ + (1 - d)Pf + a(/?/-/?g) 


or 


V 

^ m k 


Pm 


hin — 

Pm 


(A.16) 


(A. 17) 

(A.18) 
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Appendix B 


B.l Single phase perturbed enthalpy equation 

The single phase energy equation is given by (see equation 3, chapter 2) 

(B.l) 




dt ‘ dz PfA^ 

On perturbing equation (B.l)and neglecting second order terms we get 


I V I V I V' 

or sh f + Vf,, -:r- + -T~ 


‘ f ‘ ’ fss 


dz 


dz 


dz pfA, 


But at steady state from equation (B.l) 

V = 

az PfA, 

so on putting this value in equation (B.2) we get 

dh'r qlP, 


sh\ + V, .l-A + i/;— = 0 


dh\ 


or 


_£+_£_ /i'. 


/7/A Ks 


On integration of equation (B.5), we get 


h'^=e 


■J] 


"ir 






J 


-afe 


^vlsPfA 


i—Vf e liz + /f 












(B.2) 


(B.3) 


(B.4) 


(B.5) 


(B.6) 


(B.7) 
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Let Xss - steady state boiling boundary position since it is assumed that in single- 
phase region h^ = 0 


So putting at z = 0, h^ = 0 in equation (B.7), we get 


So 


so 


0 = 


A = 


-q:PhV; 

PfAcVf33S 

q:vP.,v; 

PfAeVf^^S 


+ A 




’ vl 


■q:PnV; .^7 


PfAeVf^^S 


e + 


n" P V' 

HwH 

PfAcVfe^s 


(B.8) 


(B.9) 


“ 7, aJv^s 


1-e 


B.2 Volumetric flux propagation equation 

For two phase region conservation of vapour phase 


^(2p5l+±(ap V ).r 


dz 


^apfp, ^ 


I P. 


iV- 

gJ 

m y 


(B.IO) 


for heated region 

(B.ll) 


Substituting equation (A.l 1) into equation (B.l 1) we get 


d 




® dz 

: p .| 

^ dz 


t ^ ^ g ( apjxV^Lp^ 


f Pg ’ g; 
Pm J 


dz 


(B.12) 


J 


apJJa{Pf-Pg)-pf} 


Pm 


r,-|(ap.V,) 


(B.13) 
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Since + (1 - <3r)p^ 

or Pm =P/-a(p/-pJ 

The elementary form of the mass conservation equation is obtained from equation 


(B.13)as 


or 


|(“P.)+|(«PsV.)=r. 


(B.14) 


or 


a az'' 


:r 

® dt ^ 


dz 


or ^ + ^(av): 

a az' *' 


g ^ 

Pg Pg 


^Pg 


^Pg 


dt dz 


IS 


or 

dt dz^ ® pg pg Dt 


where 


(B.15) 


D/ d^'^^^'dz 

Similarly the elementary form of mass conservation equation for the liquid phase 


or 


or 


dt 


(l - a)p^ + £ 1 = 


(B.16) 


^-|(ap,).p,|(l-aK*a-a)V,f = -r. 
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or 


or 


■Pf 


da d 
dt dz 


(l-a)v, 


= -r,-(l_a)%-(i-a)V,-?e£ 


dt dz^ Pf p 


f L 


a 


^+v,^ 

a dz 


dz 


da d 


or + A(i-a)Vf=::ii-l^^ 

dt dz p^ Dt 


Adding equations (B.17) and (B.15) together, we get 


(BAT) 


_d 

dz 


or 


((xV^,+(l-a)Vf) = r 




a Dpf ^ (1-a) Dpf 

^pg pf; 


_pg Dt Pf Dt 




PfP, 


a P>pg ^ {\-a)Dpj. 

[ 


(B.18) 


(B.19) 


^p^ Dt pj Dt 

The finst tenn on the right hand side of equation (B.19) expresses the volumetric 
source due to the phase change and the second is the volumetric sink due to 
compressibility. 

it should be noted here that if each phase undergoes isochoric process the 
compressibility effect drops. 

So, we get ^'ohx?7ici!ric flux 'propagation equation 

8j TjjA/? 


dz PfP^ 

B.3 Density Propagation Equation 

equation (B.14) can be further simplified as 


(B.20) 


(B.21) 
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da 5p d(\ 5p 
or P«‘:;r + rx--^ + p — (aV )+ 

® at a ®az^ ® az 


(av)+Va:^ = r, (B.22) 


In NCL system pressure drop is very less in comparison to operating pressure. So 
variation of Pg can be neglected (isochoric process), so 

‘^■(aVj=r, (B.23) 


a 

* a az 


+ pj 


a az. p. 


from (A. 1 5) and (A. 16) F = F . + y 




aa , aa , aj a / ) r 

or +a-— +— ^aV^j = -— 

a az az a pg 




or 


or 


da 

_ 

aof 

dt 


7 + F^ + « 


aa 


da 

dz p 


J + V^ +a- 


dV. 


^ 

da 


da 

dz 


oT^Ap 


g PfPg 

r 

aLp 




1- 


Pf 


since Ck=j + Vg 


so 


aa 

“a” 


+ C 


aa 


dz pg pf 


(B.24) 


(A.14) 

(B.25) 



_a 

' Pf ~Pm 

^ 5 

' Pf- Pm 

=;£l£=L 

or 

dt 

Ap 

-hCj. — 

dz 

Ap _ 

Pg Pf 
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or 


or 


1 

A/? 




dt 


dz 


g Pm 


r 

pgPf 


^P”' - '^tP’-^P 


dt 


On putting Q. ■ 


dz 


r.Ap 

pfpg 


PfPg 


We get ;./?/ propagation equation 


--f+o‘^- = -P.Ci 

o( oz 


For sSteady state condition 


c- A’-h. = -p q 

'"in ^ rnas‘“‘s 


for heated region 


For adiabtitic region Qss = 0-0 
so equation (B.29) can be written as 

^Pmss _ Q fQi- adiabatic region 


C 


ks% 


OZ 


B.4 Perturbed Density 

On perturbing equation (B.28) we get 


(B.26) 


(B.27) 


(B.28) 


(B.29) 


(B.30) 


|(/’™ +p:)dQ. +c;)|(p™, 

dz dz 


or -T.,, 


or sp'„ + C 


dz 


JS 3 S 


(B.31) 
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or 


On putting the value of from equation (B.29) 

(^ + n., K + Q, ^ = c; 

dz Ck„ 


or 


^Pm S + O,. , _,a^p 


p;, = Ck' 


Ck^, 


f>/. Ck,, 

Solution of equation (B.34) is 
Mixture density perturbation equation 


r 5 » lit 

f i t’l 


ik 




rc -:n. ,/r . .. f 


dz + B 


to- 


Now CK(t) - Vr, (t) + Oss (z-Mt)) + Vg,- 
since Vr, and Vgj do not change with z 
dCkss ” Qss dzi 

On substituting equation (B.36) into equation (B.35), we get 


4,1 f "L* 
p'^s:e ““ 


- 

f i+0„ dCtf 

— 

t^kPmss^ss 

J 0.-C.. 

dz + B 

J 

'-to 

- 



540, 


‘4nC 


or 


/>„ =e 


kM 


| 40 «, 


- 

j+0^ , ^ 


^S^ssPmss 

^kss 

e “ 

dz + B 


or p: =•■ (Ck,. )■ 5,: 

Now from equation (B.36) 


')S^ dz + B 

^h$ 


(B.32) 

(B.33) 

(B.34) 


(B.35) 


(B.36) 


(B.37) 
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or 


dz=^^ 

n.. 


so 


) «« 




(B.38) 




(B.29) 

or 

c, 

^kss 


or 

d/y cy 


On integration, lufi 

jck 

SS 

or 

Pnm ^ C^s ^ ^ 


or 

Ptms ( -C- constant 

(B.39) 

at 

P«|$S ^ pfs ^kss 

^hs Ms, SO from (B.39) we get 

Pnm - ■■ Pc ,xss 


or 

u 

on'-^ 

ii 

(B.40) 


Substituting equation (B.40) into equation (B.38), we get 


p>ic^ )•('§■-) 


f^dC^+B 

^kss 


pXc,. rf'S”) 




I /-'2 
'•'Jkis 


Ck 


(B.41) 
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/ 


Pm ) Pf^kss^Zss Q 


yiC+B 


(B.42) 


In equation (B.42) to find the constant B, it can be assumed that a positive perturbation in 
boiling boundary results in a negative perturbation in the mixture density. To a first 
approximation it can be written from equation (B.39) as 

= = (3 43) 


SO p, 


= ^ssPf-^ _ ^ssPf 

yTt jOt 




(B.44) 


Substitution of equation (B.44) into equation (B.42) for the condition that at 


u,..- 

PmMs ~ = (^Us ) ( ^ Pf^]!SsM 


* s-Q, 


■(^kss,;iss ) +.S 


or P. 




B = “"(’f' 




I 

s-Cl. 


(B.45) 


(B.46) 


Substituting equation (B.46) into equation (B.42), we get 


Pm = (Ckss ) ( o- ) |p/Ck^ C'k (Cto ) o„ 


+ K4 ^s,Pf^ii^)K-PA- 


P'm ~ ~Q~ P/Cks,.te Q Cte + I (Ckss..ls. )"“ [^^«P/^I s-Q 


%- (c.„ )■ 


(BAT) 
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Q 


,v - LI 


^k'is./lss 


>2 


.n V '“ks!. / 


c 


c: 




r _ ^/' 


j+Q„ 


-Cl 


kss.^s 


5-0 


pf 


V^Cj^s^jsj J 


~'kss,,iss 

V ^kss J 


since >,' ~ r, (s)Vf, 


(B.48) 


B.5 Perturbation Equation for centre of mass velocity (for core) 

On pcrturlnng both sides of equation (A.17), we get 


r - [/' -- f' , c 


PfK- 

Pmss + P'm 


^ in ^ i 


P/-^» 


Pm 


1 + 


o' ^ 

Jr m 

Pmss 


\ f^mss J 


or 


i;,. 

Pn»s 


Pr, 


\ 


1- 

V Pm^s J 


Now from equation (AJ7) for steady state condition 


or 


F 

WKj in 


E-li. 


(B.49) 


so 


K C'l + 


Pr^ffPm 


(B.50) 


PfV„. , 

or v;=r,(s)v;;+^p:„ 


(B.51) 


B.6 Perturbation of Density Propagation equation For Riser 

Density pi . .panation equation can be obtained by putting Q = 0 in equation (B.28),we 


got 


SO 


^:P« +,Q%- = o 
Pi dz 


(B.52) 
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On perturbing both sides, we get 


d 


3 

{Pmss + + C'k )— +P'„) = 0 


or spl + Cj 


'kxs 


^^ + C'te^ + Cl-^ = 0 


dz 


dz 


(B.53) 


For steady state condition equation (B.52) can be written as 


K - u 


So eqtiation (B.53) can be written 




ki% 


4- c' ^ = 0 


d= 


dz 


(B.54) 


o f 

on neglecting second order tenn C[ , we get 




(B.55) 


or 


fOm 

Pi 


S J 

——dz 


C 


ks$ 


or in plf^ “ — -}* C [Because Ckss does not change with z] 


B.C , at Z *■ L, Ckss “ ^kss.L j Pm Pm,L 
(end of channel) 


So Inpl,;. ="“ + C 


C = Inp:.,,-f' 


sL 


'hs,L 


(B.56) 


or 
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So perturbed density for riser is given by 


Pm = Pm.l 




B.7 Peturbation of centre of mass of velocity in Riser 

on perturbing equation (A 1.1 7), we get 


r 4. V' =c +c' -■ 

^ m.ft ^ ^ kss ^ 


Pm,, + P'm 


-c +r' 

Us ^ / 


PfV, 


1 + 


Pm 


or 


V t'm.sj 


r I !'' ■■■ (” I r' I Pf^sPm 


(B.57) 


(B.58) 


Pmss Pmss 

Now for steady state condition equation (A. 17) can be written as 


F rj (" 


PfK 




(B.59) 


On putting this value in equation (B.58), we get 


y -f' , Pf^gjPm 

Pmss 


(B.60) 


On putting the value of p'„ from equation (B.57) in equation (B.60), we get 


y r 4 . EiIm. n' 

~ '■ * * ”2 Pm.l. 


s(x-l) 

L W J 


pL 


(B.61) 
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B.8 Perturbation of Pressure drop components 

B.8(a) For 1-<|> unhcated regions 

(I) Gravitational component 

APk = I Pfgd:'' 

on perturbing equation (B.62), we get 
Apr.. + Ap; = I Pfgdz 


fnir steady state condition of equation (B.62), we get 

Ap,,s I 

.so Ap; 0.0 


(2) frictional component 

2D 


Ap, = I 


Vidz 


on perturbing equation (B.65), we get 

Ap..+Ap', = |^i(v,.,+v;)‘to 


■=|;--K + v;=+2v.„v;)l 

on neglecting (second order term) 

Ap,u + Ap; “ + 2Vfi„V ' )l 

Now for steady state equation (B.65) can be written as 

AP„ = f 


(B.62) 


(B.63) 


(B.64) 


(B.65) 


(B.66) 
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(B.67) 


iTorn equation (B.67) and (B.66), we get 

Ap' - I V' 

^ ^ D 

(3) Acceleration component 

Ap., 0.0 

So on perturbing equation (B.68) 

Ap... ■» Ap; == 0 

b‘or .steady state condition of (B.69) 

Ap.v. 0,0 

f rom equation (B.7()) and (B.7I), we get 

Ap; : 0.0 

(4) I Ai.ss component 

Ap, k,PfV,t’ 

or Api,, + Ap; = k,Pf + V;, )' 

or Ap(,., 4*Ap, “kjPf^Vgj, ■t"2V,5jjVg) 

For steady state condition of equation (B.73) 

Api,» ~ kjPfV,,^ 

From equation (B.74) and (B.75), we get 

Ap; « 2k,p,v,,,v; 

(5) Slip component 


(B.68) 

(B.69) 

(B.70) 

(B.71) 

(B.72) 

(B.73) 

(B.74) 

(B.75) 

(B.76) 
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Aps = 0.0 


(B.77) 


so in the same way as for perturbed acceleration component, we can get 

Ap;=0.0 (B.78) 

(6) Inertia component 

(i. dV„ 

= (B.79) 

on perturbing both sides of equation ( B.79) 

Ap,..; + Ap; - I Pr + V;;)dz 
for .steady state condition of equation (B.79) 

Ap,„ - f 


or Ap,^,^ -= 0 

so 0 + Ap[ = I PfSV^dz 

or Ap;-p,sV'L (B.80) 

B.8 (b) Perturbed components of Ap for single phase heated region (in core region) 

(I) Gravitational component 

Apg = j[pfgdz (B-81) 

On perturbing equation ( B.8 1), we get 

AP„. + Ap; = 

or Ap„4Ap; = pp,g'i2+£'“pr8'>z 
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m 


For steady state condition equation (B.82) takes the fom 

= |'“pfgdz (B.83) 


from equation {B.82) and equation (B.83), we get 


Ap; 

f“* Pfgdz 


or Apj. 



or Ap; 

-p,gX' 

(B.84) 

or A[^'. 

p,gr,(s)V' 

(B.85) 

(2) Frictional c«»mponcnt 


Ap, 

, r v=dr. 

i'l 21) * 

(B.86) 


or A/i,„ +V/ 



(B.87) 


Now frtvm equation {B.86) 


Ap, 


Elil .hence for steady state 
2D 
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On putting tfiiK value in equation (B.87), we get 

-v; “ (b.88) 

(3) Acceleration component 

Ap., 0.0 (B.89) 

So as obtained in case of single phase unheated region, we can get 

Apl 0.0 (B.90) 

(4) Loss component 

Ap, k,p,V,; (B.91) 

now in the same way as we obtained for single phase unhealed region we can 
Ap; =’k,p,V,.,VJ (B.92) 

(5) Slip C,’on 5 |x>ncnl 

Ap, 0.0 (B-93) 

So on perturbing equation (B.93), we get 

But for steady state condition of equation ( B.93) 

Ap, 0.0 

so 0 f Ap', 0,0 

St) Ap' 0.0 (B.94) 


(6) Inertia C’omponent 



Api — 0.0 ,so in the s<iiTic way as we obtained for the slip component, we can get 
Ap' -- 0.0 (B 95) 

B.9 Perturbed Pressure Drop Components in two-phase region 
B.9 (a) For the core region 
( 1 ICiravitational component 

A/’, - (B.96) 

On perturbing equation ( B.%), we get 
or Ap^,„, 4 Ap; - +p;,)gdz 

£(p ■» p;„kl/.-£“''pfgdz 

or Ap^,, + Ap; - p„«,gdz + £ p;gz - p^gdz (B.97) 

i' or steady slate of equation ( B.96) 

APp, I p„„gti/. 

On putting this value in equation (B.97) we get 

Ap; - p;,gtiz -£"'^ pfgdz (B-98) 

Ap; - £p;,gd'^'“Pfg^' 


On putting the value of p;, from equation (B.48), we get 
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dz 




(CJ-. 


■p/g^' 


dC 

On putting dz = — ^ from equation ( B.36) 

2)1 


A/?; = g 


Pf 


Qw.teQ .(Osi) ]i 


+ g 


ro. 




r ^ 


s-a 


£l 

SS J 


(a„)n„ 


za 

V ^ 






Ap; = g 


Ck„,^„Ck'-[(ck.,^)-’-(Ck„,J-'' 


S-Q 




+ A'- 


Ck' 


s-a 


PfiCkJn,, 


SS J 


o., 


V ^ 





(2) i'rictional component 

Ap, = |^P,V> 


f. 


f m „ J /-2 


ID 


P^y m '2,D ' 


On perturbing equation (B.lOl) 


Ap,„ +Ap', = £^(p™. +p:.Xv.„ +v;fdz-ip,vir 


Ap,„ + Ap', = £^(p„ +p;„XvL +2V^ +V;)i2-ip,Vr^V 

= f ^(p„v;., +2V„,p,„V; +p'„V4.)iz-ip,V4l 


(B.99) 


(B.lOO) 


(B.lOl) 


(B.102) 


' (B.103) 


80 



for steady state condition of equation (B.lOl) 


or 


Ap, - f dz 


on putting this value in equation (B.103) 

V, r 






Ilnd term 


(B.104) 


t 

D 


I.sttcnn f (p«.-. V„ . )V' dz 

J I* \r nrt^% f ^ * 




'« '-■”>• I’, v„.. I v:d/. 


(B.105) 


on putting the value of from equation (B.51) 


1st term {>, V^,,, 


f r,(5)V'd2 + f ^^tp'.dz 

*■* “ Pnas 


(B.106) 


j;pfV„j-,{s)v;;(L-xJ 

If « 

Substituting equation (B.40) into equation (B.107), we get 


1st tcmi ■- 1 2 (s)VI , (l - 4 ) + 

fii'f 

1st tenn 


A 

D 


vA 


c 


^r^Cldz 


2 

ksSyAss 


L V. 


D C 


kssMs 


- f 


(B.107) 


(B.108) 


(B.109) 
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Substituting equation ( B.48) into equation (B.109), we get 


f ' -^P^. r, (J )r; (i - 4 ) + ^ 

^ ^ ^kssMs 


K D 

a 


1 7 -S- f -^c; 




Pf 

( c \ 

c 

^kss,Ass 

^ ^ksz J 






or first tcmi - t^, V„,,r 3 (s)V'(L-Xj 


/■ F, F. 
Z™.. » /7 

V) r- 


(B.llO) 


where 


i p's'ij- 


(B.lll) 


or 


t 


Q 

.. -- o 


da 


Cl, 


Pf 

( c ^ 

'^kss,Xss 

^kssMs 

K, ^kss ) 


5+0^; 




On putting dz ^ from equation (B.36) 

tffe ^ iP 


P_Ppi>,r 


Ck' 


t 


dC 


hm 






5-0, 


1 / ^ 


kss 


a 

'^hs 


(B.112) 


F = - (Q,,, - 

5 - O.. 


U'- 


c: 


5-0 


k I, /I 

Pakss,kss 


1 

L 

la^-s 


^kss \ 


20„-5 


I 

P 

^ 1 

! 

Ass 


(B.113) 
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_ Q ^hsMs^k {^kss,L ^kssM ) 

*3 IWk M ^ « 


4 - 




c 


5 - Q, 


fP/Q«";Lsj 


/ v2n„ 5 . 

[Ckss,LT^ -(CkssMs) 


ir 


Q, 




Ilnd temi --= y-p'^V^dz 


2D 


On substituting equation (A. 17) into equation (B.l 15), we get 


f. 

lliul term 

4*' , 


U. p' 
2D 


c 


ksS 


K 


P^E 
Pms J 


dz 


/. 

2D 


i 


c 






kss 


kss 


Pf^kss,. 


P'm<i^ 


hsMs J 


/v<', = ^'P,V(i-4)+^4 


-V,„F 


kss,hs 


(B.l 14) 




(B.l 15) 


Jn 

W J*.. 


f 


/« 

2D 

L 

ID 


C V 

^kssMs J 

r, 




P'J^ 


I 




'^kss.Ms J 


Clp'Jz 






\ ^ kssMii ) 


(B.l 16) 


A 

ID 


y 


cf 


kss, Ms 


F-^^p,yD' 


(B.l 17) 
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(3) Slip component 



Pf-Pm PfPg 
,Pn, -Pg Pm 


V' 


dz 


Assuming that p,n » pg 


A/?, 


Pf-Pn 

pI 


PfP/l 


(B.118) 


On perturbing equation ( B.l 18), we get 




Pmss.L L,p 


m,L 


P 


f Pmss,L Pm,L ^ ^ T/2 


1 + 

\ PmsSyL ) 


(B.l 19) 


P f Pmss,L Pm,L 


Pm&sM 


2/^: 


m,Z. 


1- 


r mss X J 


PfP/. 


Pt ”P 




2p; 


m,L 


Pm,L 


r 


2 

r niss,L 


2 p:n,L 


M 


Ptmiuh J\ Pnt5.s.L / 

2pm,L(pf "" Pm$s,L 


1- 

Pniss,L J J 


■PfPgV, 


2 

g ' gi 


Pi P«m,L ^ ^ \r2 
PiPr^si* 




Prms.L 


) v2 Pni,LPfPgy[ 
-pfPgVg, 75 

Kmss.L 


or Ap' 


2p;a.PfPgVj ] 2p;,,LPfPgV^. Pm.LPfPgVa- 


PnwJ. 


2 

Ptos.L 


2 

r mss,L 


(B.120) 


or 


or 


. , Pm.lPfPg^lS 2pm.l.PfPgY 





pLj.pfp/a 


PmnJ, 


2p f 


PmstJ. 


(B.121) 


84 



(4) Inertia component 


Api = {l 


P m 


av„ 


a 


•dz 


On perturbing equation (B-122), we get 


A/’,,, + Ap; = 


mss ' ' 

dt 

dV 9F''^ 

\Lsb. + ' 

dt dt 


4z 


(B.122) 


(n +p' )^dz 

..r'‘ ‘ a 


or 


( since — = 0 ) 
at 


APh, 


On neglecting the second order term in (B.123), We get 

^^.sp^v;dz 


or 


AP.„ + Ap; = [ sp„„,V;dz - (“ spfV^dz 


or 


Ap.rt + Ap' = ^ sp„^,Vn,dz-spfVf,A. 


on 


neglecting second order tenn,we get 


or Ap 


,,+Ap;-|' 


(B.124) 


now 


for steady state condition 


ofequationC B.122) can be written as 
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(B.125) 


rp.»^' 32 =o 

at 

so Ap; = f sp„,,,v;ciz 


on putting the value of V; from equation (B. 51), and value of 


wc get 





its. In 


ki% 


r,(5)F; + 


'-kss 


Pf 


kssyhs 


dz 


or 



si\is)yuiz^ 


c 

v 

c, 


■P'J: 


ksSyAss 


I '*u 

now it * ’■ 

n.. 


from equation (B.36) 


so A/>, '■ j 


At/' f' 


“ ^^is^Jax ^tssjss 




,, dC, 


kss 


Q. 




Q. 


^ kss, Ms ^ss^kss,Ms 


f C„,P'JC,^ 


lcKr:jCK^l-^c.p, 


^ss r^f ^ ^kss,;iss 


c 


kss 


kss 




C 


j+n„ 


kss, Ms 

\ Ckss J 


dC, 


kss 


from equation (B.40), 
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\ ss,Jiss 




^f^ss,;us 


^ * <1.^1 )”" ^/ 


-i?i£_c;T-i^i^ f(c, ,ff‘-(c. 

s-D^ ^ 


‘M’: ^ I Pn,»sV,;,<Jz 


^^■.. ^ksU« 


' , ' "■■' , • -•-%-- In^ 

A—i-ij,, '^issMs 


c'i ^‘■-+o,.i7 — 


S-Q, 


Q,,-s 




.. PrCua.«I'2(s)K , C,,., 

Ap^ ^., In- 


F .V C' 




Pf^kssM! 


i-pAQ„A' 


£Si?» 

s ~ £ 2 „ 


r \ ( c I 

^kssMs J _][ Q 

"■ iS 4* £^*f j . ^If ^t.Xs.'! i 


(B.126) 
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(5) Acceleration tcnn 



mass conservation equation is 


'a 




(pn,Vj = 0 




flp 5V,„ 

.-'-W.v +f) — 
ft/. dz 


so 



^Pni Y 

"'a'" f>. 


(B.127) 


(B.128) 


On putting this value in equation {B.127), we get 


AP. 



^32- -V 

a ■" az J 


Oil perturbing Iwth sides of equation (B.129), we get 

^Pmss ”■ 

dt dz 


(B.129) 


Api " " .( ^ )m 



I**or steady state condition 


from equation (B.129), we get 


(B.130) 
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(B.131) 


Ap,., = ~ f' 

From equation (B.I30) and (B.131), we get 

AP', “-f V„„sp'„d2- f- Vi^^dz 

■ f -v.,„ v: -^i-'-’^idz 

K d7. 



f,. Pf^,] 

^ hi 


V / J 


[Putting equation (A. 17) for Vmss] 


(B.132) 

(B.133) 


Now on putting the value orp,a,s from equation (B.40) into equation (B.133), we get 


I . I 


v.. 




sp;dz 




f ’ 


I Ci„p:d-/ 


I 


c\ 


4 / 


4 n„A' 


I w 


-r 




c 

TI o c ^ 




iifk" 

1 Pf 


J 

l J 



dz 


on putting value of dz from equation (B.36) 


l- 


('k 




Q„ ^ 
s-'Q. 


<PfC, 


t. 


dC, 


kss 


K Q„C 


kss 
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+ 


f y \ 

V ^ < n. «i / "“t.t 


-.1 /- cfC 


n /I' 
a ' “ 


s-Q, 




less 


ic.J- 


r., “1 

S 

c, 

1 — '•' . 

t ' 

\ A / 

.. " fC~ 

In 

. ^kss^ss_ 


/ ^ 




A \'s,.iu / 


c 


a’-Q„ 


W(Q- 




n„ 






('\..., i't' (t\„..J“^ 


';■ ' 


! ( 


I . f n / 


.V £2. 




ISK 


ln-^‘-i!i 

^ * ss,Xss 


(' 

,V £2. 


- \PrKJi 


^SS_ 

Pss-^^ 




mm 


1 


V' . &/. 

”’■■ (T/, 


D* V^ do' 

r^mn ^ 


Prms^niss 


“ pL ^ 


(B.134) 

(B.135) 


now V,„v( pi Vf„x from conservation of mass and putting the value of 


Ml 

dz 


equation (B.34), we get 



( 


&’'k 

i 

Pm%\ 


A 4- £2, 


C 




) 


P^ssPm 

~c/ 


■CJz 


ks$ 


on putting the value of pn^-i from equation (B.40), we get 






dz 


P f^kssMs ^kss 


from 
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(on putting the value of from eqn. (B.48)) 


r' j ( Ci c \ f C 

'^kssMs r" r r r> T ^ss'-'k _ ac 




^“c«- ^Jtcc JLw / \ SS J 


Ms 
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(B.138) 


+ a„A'- 


Q 


5 — Q 


« f'f 

'^k 

ss J 


v^r Q, 




Cl,, -s 


kss,L 


-{Ck 


0.gs \^kss^^s 



+ 


pn 


J Pf^kss.^ss 


“ CM 


^ c ^ 

^kss,L 
\^^kssMs ) 


^ fiss 


s + Cl,, 

C kssMs ) 




c, 


kss,L 


-Cl 


C 


kssM^ 


+ 


^r^ Q' —^hl—C 
S-Cl ' 


-"ss / 




PfCfiss^^ ^to,Z^ 
Xss P'kss,kss 


now f 2 V„^ 5 V„ 

•^ss 




dz 


< 


'^Pf^fiss 


CM^^P'm 

Pmss 


Pmss^ss_ ^ 2 : 

a 


'kss 


(Putting equation (B.50) and equation (B.29)) 


Pf^.- r 1 

-2pA^^^ssLlcM-y-p^ ^ 


Pn 


kss 


(on putting equation (B.36) for dz) 


f ^^/iss^g/ f n' dC 

: -lp,Vf,.C I 


kss 


fiss^k , 

^ kss Ms kss Ms 


(B.139) 


(B.140) 


(B.141) 


(B.142) 
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f c^,pyc, 


+ (nr- 


5-Q, 


c; W(c, 


kss^Zss y " 


0„-5 


(c.„.)%r-(c,„ 


Ms) 


(B.137) 


so f 2)>;„ VI dz = -2y0 C[ In 

” az * a„ a! 


(B.143) 


(6) Loss component 
Ap,=K,p„V^ 


(B.144) 


aa „ ^^^p',=KXp„,,+p'Jr^ss^rJ =K,{p^.+p'Jyl+2v„ji) 

= K,p^, + 2 K,p^, v; 4- K.p;v^, 

From equation (B.144) for steady state condition 


Apus =K,P^ 3 V^ 


(B.145) 


so i^\-lK,p^sV,^rm^K,p'Jl 


(B.146) 


on putting the value of V^mi,we get 


^K.p^r^v: +K,p\ 


tip fissV fiss 


Now, on putting the value of p^^ from equation (B.40), we get 




(B.147) 
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B.9 (b) Perturbed pressure drop for riser (unheated 2-(j> region) 

(1) Gravitational component 

^Pg = fpmgdz (B.148) 

on perturbing both sides of equation (B.148) we get 

^Pgss+Ap; = f (p„^3+Pn,^dz 

or + Ap'g = {" Pn^s gdz + jf" p'„ gdz (B. 149) 

for steady state condition equation (B.148) can be written as 

= I' P^ssSdz (B.150) 

from equation (B.149) and equation (B.150), we get 

APj = gdz 

on putting the value of from equation (B.57), we get 


= g 



s{z-L) 

e dz 


^P'g = P'n.,Lg 


c ^ 


Cks!.L 


- 1.0 


(B.151) 


(2) Frictional component 

AP,= {"^P„V> (B.152) 

perturbing both sides, we get 
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Ap,„ +Ap', = f ^(p™ +p'Jv.„ +v;)=dz 
= +p'.Xvl + 2v^v;)te 


Ap,„ +Ap; = |■'iv4.dz^-f2p^iv^v:d2 


+ (^'p' iis-V^ dz 

i ?n 


2D 

for steady state condition equation (B.152) can be written as 
'L 


(B.153) 


APrs. = f ^JPn..V^3dz 


(B.154) 


from equation (B.153) and (B.154), we get 

Ap; = |"§-(p™,V„)v:d2+fip'„Vi.dz 


= ini.(p V )f''V'dz+— f"’ d' dz 


or on putting value of andp'„, from equation (B.61) and (B.57) respectively, we get 


Ap} 


s{z-L) 


r^t , P f^sPm,L ^ 
^ifc 2 


\dz 


j-£sLV^ 
D ^ 




s(,z-L) 


dz 


or Ap',=^{p,„r,JC,(L,-L) + 


Pf^sPm,L 


7 (7 ^ 

V S j 


PI 
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(B.155) 




• 1.0 


+ 


JjS-Y^ n' 

2Z) ^ 


( 


V ^ J 


(1) Slip component 

Ap. = f 


^ Pf Pm PfPg y2 
5z[pm~Pg Pm 


!dz 


or 


APs = 


Pf Pm PfPg y2 
P m P g Pm 


On perturbing both sides, we get 
^Ps + ^p1 = 


iL, 


Pf Pmss Pni p. ^ 

PfPg ^gj 


(p mss P m ) 


since pm » pg 


pf Pmss P 

,(pL + 2p^5p'„,) 


TlPfPgV^ 


P f P mss P n 


2 

P mss 


i-^)pfp,v= 

V Pmss ) 


Pf P mss 

2 

J 2p „ 

PfPgVg- 

Pmss 

\ P mss > 





nL, 


PfP.V, 


2 

g ’ gJ 


JL 


Pmss \ Pmss / 

now for steady state condition of equation (B.156) 


Ap, 


Pf Pniss PfPg y2 
_ ''a 

L P mss P g P mss J i 




'A-Q 1 
- 1.0 1 


(B.156) 


(B.157) 


(B.158) 
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On neglecting 2"^ order term 


Ap; = J 

H ntss 0 ^ 

(B.163) 

Pmss^ ^ 

t'v:dz 



on putting the value of from equation (B.61) 




c + 


Pf^Spm,L^ 


^kss L 


P mss, L 


\dz 


^p'i = 


c;(.,)(A -L)+ 




2 

r mss, L 


C, 


kss,L 


sm-L) 

^ks5,L 


- 1.0 


5. Acceleration component 

Ap„ = f 


now from conservation of mass equation 


gpm 

a 


+— (p V ) = o 

oz 


or 


a ^ & 


on putting equation (B.166) into equation (B.165) 



(B.164) 


(B.165) 


(B.166) 


or 


dz 


(B.167) 



on perturbing equation (B.167), we get 


Ap„ +Ap’. =-{"(v„ +v;)|.(p^ + P ',) cl2 


dt 


-fc 


oz 


(B.168) 


Ap..+Ap'. =-f‘V^sp;d2-['Vi,:??™d2 


dz 


|-y 2 5 p; 


w 5p 


dz-{"2V^V;A--d2 


dz ^ az 

now for steady state condition of equation (B.52) 


^Pmss A 

'-'te “ ^ 


dz 


(B.169) 


or 


ap 


mss __ 


az 


0 


Also for steady state condition of equation (B.167), we get 

'-dz 




dz 


(B.170) 


(B.171) 


on putting equation (B.170) and equation (B.171) into equation (B.169), we get 


Ap’. = - f s p;dz - f ■^<*2 


-t 


sU~l) 

V s p' ,e — f* 

^ mss 


VL — dz 

OZ 




^ c ^ 

^kssX 


sm-i) 


- 1.0 
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Icfz 


c 


kss J 


^mss ^ Pm 


Z' (7 ^ 

^lcss,L 


s(Li-L) 

e ~1.0 


+ 


s 

mss 


a 


kss 


I' Pm^-Z 


V o' C 

^ mss r' mj. ks% X 


s{h-l) 

\ x r. 


- 1.0 


+ ■ 


s 

^ mss ^ 


s(z-L) 


Ck, 


ssX 


fv> 


dz 


^mss PmX ^kssX f 


s(Xrn 

^kis,L 


- 1.0 




o' 

^ mssr mX 


C 


kssX 


c ^ 

'^hs,L 



6. Loss term 


Ap,=K,p„V^ 

now in the same way as we obtained for core region, we can get 

Apj = ^ “^K.jpj^^n^mss 


^kss L 10 


(B.172) 


(B.173) 


(B.174) 
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